Visualizing cover data at each step of the data preparation workflow

library(tidyverse) 
library(sf) 
library(terra)
library(ggpubr)
ggplot2::theme_set(theme_classic())

This document visualizes the data of vegetation cover by functional type that we use in subsequent analysis at each step of the data-processing workflow. These maps serve both as a reference and a check that the workflow is performing as intended.

For the sake of simplified comparison, I show data from 2016, 2020, and 2022, which are years for which we have easy to access rasters of FIA data to compare to.

RAP

Initial dataset

## Reading layer `vegCompPoints' from data source 
##   `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_processed/CoverData/DataForAnalysisPoints' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1266274 features and 30 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -124.7004 ymin: 24.68798 xmax: -67.10236 ymax: 49.3489
## Geodetic CRS:  WGS 84

After removing burned area

# Read in data 
dat_2 <- readRDS(file = "../../../Data_processed/CoverData/dataForAnalysis_fireRemoved.rds")
# trim to be only later than 2000, which is what we want to use for further analysis
dat_2 <- dat_2 %>% 
  filter(Year >= 2000)

## make figure for RAP 
# get data just for RAP 
RAP_2 <- dat_2 %>% 
  filter(Source == "RAP") %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))
 
ggarrange(ggplot(RAP_2 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas - RAP"),
  
  ggplot(RAP_2) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas, \nall years - RAP "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After ‘coarsifying’ LANDFIRE data

# Read in data 
dat_3 <- readRDS(file = "../../../Data_processed/CoverData/data_beforeSpatialAveraging_sampledLANDFIRE.rds")
# trim to be only later than 2000, which is what we want to use for further analysis
dat_3 <- dat_3 %>% 
  filter(Year >= 2000)

## make figure for RAP 
# get data just for RAP 
RAP_3 <- dat_3 %>% 
  filter(Source == "RAP") %>% 
  pivot_longer(cols = c(#ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC
    ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))
 
ggarrange(ggplot(RAP_3 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data - RAP"),
  
  ggplot(RAP_3) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data, \nall years - RAP "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

AIM (BLM Assessment, Inventory and Monitoring data)

Initial Dataset

## make figure for LDC 
# get data just for LDC 
LDC_1 <- dat_1 %>% 
  st_zm() %>% 
  filter(Source == "LDC") %>% 
  #st_drop_geometry() %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues)) 
  #filter(coverType == "ShrbCvr" & Year == 2016) %>% 

ggarrange(ggplot(LDC_1 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data - AIM"),
  
  ggplot(LDC_1) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data, \nall years - AIM "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After removing burned area

## make figure for LDC 
# get data just for LDC 
AIM_2 <- dat_2 %>% 
  filter(Source == "LDC") %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))

ggarrange(ggplot(AIM_2 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas - AIM"),
  
  ggplot(AIM_2) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas, \nall years - AIM "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After ‘coarsifying’ LANDFIRE data

## make figure for AIM 
# get data just for AIM 
AIM_3 <- dat_3 %>% 
  filter(Source == "LDC") %>% 
  pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))
 
ggarrange(ggplot(AIM_3 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data - AIM"),
  
  ggplot(AIM_3) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data, \nall years - AIM "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

FIA (USFS Forest Inventory and Analysis data)

Initial Dataset

## make figure for FIA 
# get data just for FIA 
 FIA_1 <- dat_1 %>% 
  st_zm() %>% 
  filter(Source == "FIA") %>% 
  #st_drop_geometry() %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues)) 

ggarrange(ggplot(FIA_1 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data - FIA"),
  
  ggplot(FIA_1) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data, \nall years - FIA "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After removing burned area

## make figure for FIA 
# get data just for FIA 
FIA_2 <- dat_2 %>% 
  filter(Source == "FIA") %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))

ggarrange(ggplot(FIA_2 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas - FIA"),
  
  ggplot(FIA_2) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas, \nall years - FIA "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After ‘coarsifying’ LANDFIRE data

## make figure for FIA 
# get data just for FIA 
FIA_3 <- dat_3 %>% 
  filter(Source == "FIA") %>% 
  pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))
 
ggarrange(ggplot(FIA_3 %>% 
  filter(Year %in% c(2016, 2020, 2022))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data - FIA"),
  
  ggplot(FIA_3) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data, \nall years - FIA "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

LANDFIRE Reference database (use different years for LANDFIRE due to data availability)

Initial Dataset

## make figure for LANDFIRE 
# get data just for LANDFIRE 
 LANDFIRE_1 <- dat_1 %>% 
  st_zm() %>% 
  filter(Source == "LANDFIRE") %>% 
  #st_drop_geometry() %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues)) 

 ggarrange(ggplot(LANDFIRE_1 %>% 
  filter(Year %in% c(2003, 2007, 2015))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data - LANDFIRE"),
  
  ggplot(LANDFIRE_1) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Initial Data, \nall years - LANDFIRE "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After removing burned area

## make figure for LANDFIRE 
# get data just for LANDFIRE 
LANDFIRE_2 <- dat_2 %>% 
  filter(Source == "LANDFIRE") %>% 
  pivot_longer(cols = c(ShrbCvr, TtlTrCv, TtlHrbC, BrGrndC#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues)) 

ggarrange(ggplot(LANDFIRE_2 %>% 
  filter(Year %in% c(2003, 2007, 2015))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas - LANDFIRE"),
  
  ggplot(LANDFIRE_2) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("Filtered Burned Areas, \nall years - LANDFIRE "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After ‘coarsifying’ LANDFIRE data

## make figure for LANDFIRE 
# get data just for LANDFIRE 
LANDFIRE_3 <- dat_3 %>% 
  filter(Source == "LANDFIRE") %>% 
  pivot_longer(cols = c(ShrubCover, TotalTreeCover, TotalHerbaceousCover, BareGroundCover#, AngTrC_, CnfTrC_, C3GrmC_, C4GrmC_, FrbCvr_
                        ), 
               names_to = "coverType", 
               values_to = "coverValues"
                ) %>% 
  filter(!is.na(coverValues))
 
ggarrange(ggplot(LANDFIRE_3 %>% 
  filter(Year %in% c(2003, 2007, 2015))) + 
  facet_grid(rows = vars(coverType), cols = vars(Year)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data - LANDFIRE"),
  
  ggplot(LANDFIRE_3) + 
  facet_grid(rows = vars(coverType)) + 
  stat_summary_2d(aes(x = Lon, y = Lat, z = coverValues), fun = mean, binwidth = .1) + 
     scale_fill_viridis_c(option = "A", guide = guide_colorbar(title = "% cover"), 
                          limits = c(0,100)) +
    #geom_raster(aes(x = Lon, y = Lat, fill = coverValues)) + 
    #geom_point(aes(x = Lon, y = Lat, col = coverValues)) +
  xlim(c(-130,-65)) + 
    ylim(c(25, 50)) + 
  ggtitle("After 'coarsifying' LANDFIRE data, \nall years - LANDFIRE "),
  
  ncol = 2, align = "h",
  common.legend = TRUE,
  widths = c(.7, .3)
)

After data has been spatially averaged (can no longer look at data by data source)

After initial spatial averaging

## After adding climate/weather data ## After adding ecoregion data

After adding soils data and removing RAP points that are in ag/developed areas